Fuente: https://github.com/rafael-fuente/Ideal-Gas-Simulation-To-Verify-Maxwell-Boltzmann-distribution
import numpy as np
class Particle:
"""Define physics of elastic collision."""
def __init__(self, mass, radius, position, velocity):
"""Initialize a Particle object
mass the mass of particle
radius the radius of particle
position the position vector of particle
velocity the velocity vector of particle
"""
self.mass = mass
self.radius = radius
# last position and velocity
self.position = np.array(position)
self.velocity = np.array(velocity)
# all position and velocities recorded during the simulation
self.solpos = [np.copy(self.position)]
self.solvel = [np.copy(self.velocity)]
self.solvel_mag = [np.linalg.norm(np.copy(self.velocity))]
def compute_step(self, step):
"""Compute position of next step."""
self.position += step * self.velocity
self.solpos.append(np.copy(self.position))
self.solvel.append(np.copy(self.velocity))
self.solvel_mag.append(np.linalg.norm(np.copy(self.velocity)))
def check_coll(self, particle):
"""Check if there is a collision with another particle."""
r1, r2 = self.radius, particle.radius
x1, x2 = self.position, particle.position
di = x2-x1
norm = np.linalg.norm(di)
if norm-(r1+r2)*1.4 < 0:
return True
else:
return False
def compute_coll(self, particle, step):
"""Compute velocity after collision with another particle."""
m1, m2 = self.mass, particle.mass
r1, r2 = self.radius, particle.radius
v1, v2 = self.velocity, particle.velocity
x1, x2 = self.position, particle.position
di = x2-x1
norm = np.linalg.norm(di)
if norm-(r1+r2)*1.1 < step*abs(np.dot(v1-v2, di))/norm:
self.velocity = v1 - 2. * m2/(m1+m2) * np.dot(v1-v2, di) / (np.linalg.norm(di)**2.) * di
particle.velocity = v2 - 2. * m1/(m2+m1) * np.dot(v2-v1, (-di)) / (np.linalg.norm(di)**2.) * (-di)
def compute_refl(self, step, size):
"""Compute velocity after hitting an edge.
step the computation step
size the medium size
"""
r, v, x = self.radius, self.velocity, self.position
projx = step*abs(np.dot(v,np.array([1.,0.])))
projy = step*abs(np.dot(v,np.array([0.,1.])))
if abs(x[0])-r < projx or abs(size-x[0])-r < projx:
self.velocity[0] *= -1
if abs(x[1])-r < projy or abs(size-x[1])-r < projy:
self.velocity[1] *= -1.
def solve_step(particle_list, step, size):
"""Solve a step for every particle."""
# Detect edge-hitting and collision of every particle
for i in range(len(particle_list)):
particle_list[i].compute_refl(step,size)
for j in range(i+1,len(particle_list)):
particle_list[i].compute_coll(particle_list[j],step)
# Compute position of every particle
for particle in particle_list:
particle.compute_step(step)
################################################################################################################################
def init_list_random(N, radius, mass, boxsize):
"""Generate N Particle objects in a random way in a list."""
particle_list = []
for i in range(N):
# v_mag = np.random.rand(1)*6 # CAMBIE ESTO
v_mag = 3
# v_ang = np.random.rand(1)*2*np.pi # CAMBIE ESTO
v_ang = np.random.rand(1)*0.5*np.pi + np.pi
v = np.append(v_mag*np.cos(v_ang), v_mag*np.sin(v_ang))
collision = True
while(collision == True):
collision = False
# pos = radius + np.random.rand(2)*(boxsize-2*radius)
pos = radius + np.random.rand(2)*(0.5*boxsize-2/2*radius) ### CAMBIE ESTA
newparticle = Particle(mass, radius, pos, v)
for j in range(len(particle_list)):
collision = newparticle.check_coll( particle_list[j] )
if collision == True:
break
particle_list.append(newparticle)
return particle_list
particle_number = 300
boxsize = 200.
# You need a larger tfin and stepnumber to get the equilibrium state. But the computation takes more time.
tfin = 130#*30
stepnumber = 2000#*30
timestep = tfin/stepnumber
particle_list = init_list_random(particle_number, radius = 1.5, mass = 1, boxsize = 200)
# Compute simulation (It takes some time if stepnumber and particle_number are large)
for i in range(stepnumber):
print(i)
solve_step(particle_list, timestep, boxsize)
#print(i)
################################################################################################################################
# Visualization of the solution with matplotlib. It use a slider to change the time
# import matplotlib.pyplot as plt
# fig = plt.figure(figsize=(12, 6))
# ax = fig.add_subplot(1,2,1)
# hist = fig.add_subplot(1,2,2)
# plt.subplots_adjust(bottom=0.2,left=0.15)
# ax.axis('equal')
# ax.axis([-1, 30, -1, 30])
# ax.xaxis.set_visible(False)
# ax.yaxis.set_visible(False)
# ax.set_xlim([0,boxsize])
# ax.set_ylim([0,boxsize])
# # Draw Particles as circles
# circle = [None]*particle_number
# for i in range(particle_number):
# circle[i] = plt.Circle((particle_list[i].solpos[0][0], particle_list[i].solpos[0][1]), particle_list[i].radius, ec="black", lw=1.5, zorder=20)
# ax.add_patch(circle[i])
# # Graph Particles speed histogram
# vel_mod = [particle_list[i].solvel_mag[0] for i in range(len(particle_list))]
# hist.hist(vel_mod, bins= 30, density = True, label = "Simulation Data")
# hist.set_xlabel("Speed")
# hist.set_ylabel("Frecuency Density")
# Compute 2d Boltzmann distribution
#total energy should be constant for any time index
def total_Energy(particle_list, index):
return sum([particle_list[i].mass / 2. * particle_list[i].solvel_mag[index]**2 for i in range(len(particle_list))])
E = total_Energy(particle_list, 0)
Average_E = E/len(particle_list)
k = 1.38064852e-23
T = 2*Average_E/(2*k)
m = particle_list[0].mass
v = np.linspace(0,10,120)
fv = m*np.exp(-m*v**2/(2*T*k))/(2*np.pi*T*k)*2*np.pi*v
# hist.plot(v,fv, label = "Maxwell–Boltzmann distribution")
# hist.legend(loc ="upper right")
# from matplotlib.widgets import Slider
# slider_ax = plt.axes([0.1, 0.05, 0.8, 0.05])
# slider = Slider(slider_ax, # the axes object containing the slider
# 't', # the name of the slider parameter
# 0, # minimal value of the parameter
# tfin, # maximal value of the parameter
# valinit=0, # initial value of the parameter
# color = '#5c05ff'
# )
# def update(time):
# i = int(np.rint(time/timestep))
# #ax.set_title('Energy =' + str(Energy[i]))
# # Draw Particles as circles
# for j in range(particle_number):
# circle[j].center = particle_list[j].solpos[i][0], particle_list[j].solpos[i][1]
# hist.clear()
# # Graph Particles speed histogram
# vel_mod = [particle_list[j].solvel_mag[i] for j in range(len(particle_list))]
# hist.hist(vel_mod, bins= 30, density = True, label = "Simulation Data")
# hist.set_xlabel("Speed")
# hist.set_ylabel("Frecuency Density")
# # Compute 2d Boltzmann distribution
# E = total_Energy(particle_list, i)
# Average_E = E/len(particle_list)
# k = 1.38064852e-23
# T = 2*Average_E/(2*k)
# m = particle_list[0].mass
# v = np.linspace(0,10,120)
# fv = m*np.exp(-m*v**2/(2*T*k))/(2*np.pi*T*k)*2*np.pi*v
# hist.plot(v,fv, label = "Maxwell–Boltzmann distribution")
# hist.legend(loc ="upper right")
# slider.on_changed(update)
# plt.show()
pos_x_list = [particle_list[j].solpos[0][0] for j in range(len(particle_list))]
pos_x_list = np.array(pos_x_list)
pos_x_list.shape
pos_x_hist = np.histogram(pos_x_list, bins=10, range=None)
pos_x_hist
entropy_x = np.sum( pos_x_hist[0]*np.log(pos_x_hist[0]) )
entropy_x
pos_y_list = [particle_list[j].solpos[0][1] for j in range(len(particle_list))]
pos_y_hist = np.histogram(pos_y_list, bins=10, range=None)
pos_y_hist
entropy_y = np.sum( pos_y_hist[0]*np.log(pos_y_hist[0]) )
entropy_y
particle_list[0].solpos[0:5]
pos_y_hist[0]
entropies_y = []
for k in range(len(particle_list[0].solpos)):
pos_y_list = [particle_list[j].solpos[k][1] for j in range(len(particle_list))]
pos_y_hist = np.histogram(pos_y_list, bins=50, range=(0,200))
entropy_y = -np.sum( pos_y_hist[0]*np.log(pos_y_hist[0]+1) )
# entropy_y = -np.sum( np.log(pos_y_hist[0]**pos_y_hist[0]) )
entropies_y.append(entropy_y)
entropies_y = np.array(entropies_y)
entropies_y[0:5]
# print(entropy_y)
plt.figure(figsize=(12,7))
plt.plot(entropies_y)
entropies_x = []
for k in range(len(particle_list[0].solpos)):
pos_x_list = [particle_list[j].solpos[k][0] for j in range(len(particle_list))]
pos_x_hist = np.histogram(pos_x_list, bins=50, range=(0,200))
entropy_x = -np.sum( pos_x_hist[0]*np.log(pos_x_hist[0]+1) )
# entropy_x = -np.sum( np.log(pos_x_hist[0]**pos_x_hist[0]) )
entropies_x.append(entropy_x)
entropies_x = np.array(entropies_x)
entropies_x[0:5]
# print(entropy_y)
# plt.figure(figsize=(12,7))
# plt.plot(entropies_x)
# plt.figure(figsize=(12,7))
# plt.plot(entropies_x, entropies_y)
plt.figure(figsize=(12,7))
plt.plot(entropies_x + entropies_y)
vel_x_list = [particle_list[j].solvel[15][0] for j in range(len(particle_list))]
vel_x_list[0:5]
vel_x_hist = np.histogram(vel_x_list, bins=10, range=(0,10))
vel_x_hist
entropies_vx = []
for k in range(len(particle_list[0].solpos)):
vel_x_list = [particle_list[j].solvel[k][0] for j in range(len(particle_list))]
vel_x_hist = np.histogram(vel_x_list, bins=50, range=(-5,5))
entropy_vx = -np.sum( vel_x_hist[0]*np.log(vel_x_hist[0]+0.01) )
entropies_vx.append(entropy_vx)
# vel_mod = [particle_list[j].solvel_mag[i] for j in range(len(particle_list))]
# vel_hist = np.histogram(vel_mod, bins=10, range=None, normed=None, weights=None, density=True)
entropies_vx = np.array(entropies_vx)
entropies_vx[0:5]
# print(entropy_y)
entropies_vy = []
for k in range(len(particle_list[0].solpos)):
vel_y_list = [particle_list[j].solvel[k][1] for j in range(len(particle_list))]
vel_y_hist = np.histogram(vel_y_list, bins=50, range=(-5,5))
entropy_vy = -np.sum( vel_y_hist[0]*np.log(vel_y_hist[0]+0.01) )
entropies_vy.append(entropy_vy)
# vel_mod = [particle_list[j].solvel_mag[i] for j in range(len(particle_list))]
# vel_hist = np.histogram(vel_mod, bins=10, range=None, normed=None, weights=None, density=True)
entropies_vy = np.array(entropies_vy)
entropies_vy[0:5]
# print(entropy_y)
entropies_vy
np.histogram(vel_y_list, bins=10)
Entropia_r = entropies_x + entropies_y #+ entropies_vx + entropies_vy
Entropia_v = entropies_vx + entropies_vy
plt.figure(figsize=(12,7))
plt.plot(Entropia_r)
plt.show()
plt.figure(figsize=(12,7))
plt.plot(Entropia_v)
plt.show()
Entropia = Entropia_r + Entropia_v
plt.figure(figsize=(12,7))
plt.plot(Entropia)
plt.show()
# for k in range(stepnumber):
# for k in range(0, stepnumber, 50):
k = 1#stepnumber
fig = plt.figure(figsize=(12, 12))
ax = fig.add_subplot(2,2,1)
# hist = fig.add_subplot(1,2,2)
plt.subplots_adjust(bottom=0.2,left=0.15)
ax.axis('equal')
ax.axis([-1, 30, -1, 30])
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
ax.set_xlim([0,boxsize])
ax.set_ylim([0,boxsize])
# Draw Particles as circles
circle = [None]*particle_number
for i in range(particle_number):
circle[i] = plt.Circle((particle_list[i].solpos[k][0], particle_list[i].solpos[k][1]), particle_list[i].radius, ec="black", lw=1.5, zorder=20)
ax.add_patch(circle[i])
vel_y_list = [particle_list[j].solvel[k][1] for j in range(len(particle_list))]
vel_x_list = [particle_list[j].solvel[k][0] for j in range(len(particle_list))]
pos_y_list = [particle_list[j].solpos[k][1] for j in range(len(particle_list))]
pos_x_list = [particle_list[j].solpos[k][0] for j in range(len(particle_list))]
pos_partic_j = np.array(particle_list[31].solpos)
plt.plot(pos_partic_j.T[0], pos_partic_j.T[1])
pos_partic_j = np.array(particle_list[75].solpos)
plt.plot(pos_partic_j.T[0], pos_partic_j.T[1])
pos_partic_j = np.array(particle_list[42].solpos)
plt.plot(pos_partic_j.T[0], pos_partic_j.T[1])
ax.quiver(pos_x_list, pos_y_list, vel_x_list, vel_y_list, color='orange')
hist = fig.add_subplot(2,2,2)
vel_mod = [particle_list[j].solvel_mag[k] for j in range(len(particle_list))]
hist.hist(vel_mod, bins= 30, density = True, label = "Simulation Data")
hist.set_xlabel("Speed")
hist.set_ylabel("Frecuency Density")
# Compute 2d Boltzmann distribution
E = total_Energy(particle_list, k)
Average_E = E/len(particle_list)
kb = 1.38064852e-23
T = 2*Average_E/(2*kb)
m = particle_list[0].mass
v = np.linspace(0,10,120)
fv = m*np.exp(-m*v**2/(2*T*kb))/(2*np.pi*T*kb)*2*np.pi*v
hist.plot(v,fv, label = "Maxwell–Boltzmann distribution")
hist.legend(loc ="upper right")
hist_x = fig.add_subplot(2,2,3)
pos_x_list = [particle_list[j].solpos[k][0] for j in range(len(particle_list))]
hist_x.hist(pos_x_list, bins= 50, range=(0,200), density = True, label = "Simulation Data")
hist_y = fig.add_subplot(2,2,4)
pos_y_list = [particle_list[j].solpos[k][1] for j in range(len(particle_list))]
hist_y.hist(pos_y_list, bins= 50, range=(0,200), density = True, label = "Simulation Data")
# plt.quiver(x_pos, y_pos, final_vector[0], final_vector[1], scale=1e9, units="xy")
fig, ax = plt.subplots(figsize=(12, 12))
for k in range(0, stepnumber, 5): #stepnumber
ax.clear()
plt.clf()
# plotting line
# ax.plot(x, np.sinc(x**2) + np.sin(x + 2 * np.pi / duration * t), lw = 3)
# ax.set_ylim(-1.5, 2.5)
# returning numpy image
# return mplfig_to_npimage(fig)
# hist = fig.add_subplot(1,2,2)
# plt.subplots_adjust(bottom=0.2,left=0.15)
ax = fig.add_subplot(2,2,1)
# hist = fig.add_subplot(1,2,2)
plt.subplots_adjust(bottom=0.2,left=0.15)
ax.axis('equal')
ax.axis([-1, 30, -1, 30])
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
ax.set_xlim([0,boxsize])
ax.set_ylim([0,boxsize])
# Draw Particles as circles
circle = [None]*particle_number
for i in range(particle_number):
circle[i] = plt.Circle((particle_list[i].solpos[k][0], particle_list[i].solpos[k][1]), particle_list[i].radius, ec="black", lw=1.5, zorder=20)
ax.add_patch(circle[i])
vel_y_list = [particle_list[j].solvel[k][1] for j in range(len(particle_list))]
vel_x_list = [particle_list[j].solvel[k][0] for j in range(len(particle_list))]
pos_y_list = [particle_list[j].solpos[k][1] for j in range(len(particle_list))]
pos_x_list = [particle_list[j].solpos[k][0] for j in range(len(particle_list))]
pos_partic_j = np.array(particle_list[31].solpos)
plt.plot(pos_partic_j.T[0][:k], pos_partic_j.T[1][:k])
pos_partic_j = np.array(particle_list[75].solpos)
plt.plot(pos_partic_j.T[0][:k], pos_partic_j.T[1][:k])
pos_partic_j = np.array(particle_list[42].solpos)
plt.plot(pos_partic_j.T[0][:k], pos_partic_j.T[1][:k])
ax.quiver(pos_x_list, pos_y_list, vel_x_list, vel_y_list, color='orange')
hist = fig.add_subplot(2,2,2)
vel_mod = [particle_list[j].solvel_mag[k] for j in range(len(particle_list))]
hist.hist(vel_mod, bins= 12, density = True, label = "Simulation Data")
hist.set_xlabel("Speed")
hist.set_ylabel("Frecuency Density")
# Compute 2d Boltzmann distribution
E = total_Energy(particle_list, k)
Average_E = E/len(particle_list)
kb = 1.38064852e-23
T = 2*Average_E/(2*kb)
m = particle_list[0].mass
v = np.linspace(0,10,120)
fv = m*np.exp(-m*v**2/(2*T*kb))/(2*np.pi*T*kb)*2*np.pi*v
hist.plot(v,fv, label = "Maxwell–Boltzmann distribution")
hist.legend(loc ="upper right")
plt.ylim([0, 0.4])
hist_x = fig.add_subplot(2,2,3)
pos_x_list = [particle_list[j].solpos[k][0] for j in range(len(particle_list))]
hist_x.hist(pos_x_list, bins= 15, range=(0,200), density = True, label = "Simulation Data")
plt.title('Histograma de x')
plt.ylim([0, 0.02])
# hist_y = fig.add_subplot(2,2,4)
# pos_y_list = [particle_list[j].solpos[k][1] for j in range(len(particle_list))]
# hist_y.hist(pos_y_list, bins= 20, range=(0,200), density = True, label = "Simulation Data")
# plt.ylim([0, 0.02])
entrop_plot = fig.add_subplot(2,2,4)
entrop_plot.plot(Entropia[:k])
plt.title('Entropía')
plt.ylim([Entropia.min(), Entropia.max()])
plt.xlim([0, stepnumber])
plt.xlabel('Pasos de tiempo')
plt.savefig(f'gas_t_{k}.png', dpi=100)
!pip install moviepy
from moviepy.editor import ImageSequenceClip
clip = ImageSequenceClip([f'gas_t_{k}.png' for k in range(0, stepnumber, 5)], fps=20) #stepnumber
clip.ipython_display(fps = 20, loop = True, width=720)